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A HIGHER-ORDER GRADIENT FLOW SCHEME FOR A SINGULAR 
ONE-DIMENSIONAL DIFFUSION EQUATION 

BERTRAM DURING, PHILIPP FUCHS, AND ANSGAR JUNGEL 


Abstract. A nonlinear diffusion equation, interpreted as a Wasserstein gradient flow, 
is numerically solved in one space dimension using a higher-order minimizing movement 
scheme based on the BDF (backward differentiation formula) discretization. In each time 
step, the approximation is obtained as the solution of a constrained quadratic minimization 
problem on a finite-dimensional space consisting of piecewise quadratic basis functions. 
The numerical scheme conserves the mass and dissipates the G-norm of the two-step BDF 
time approximation. Numerically, also the discrete entropy and variance are decaying. 
The decay turns out to be exponential in all cases. The corresponding decay rates are 
computed numerically for various grid numbers. 


1. Introduction 

The aim of this paper is to propose and study a fully discrete higher-order variant of the 
minimizing movement scheme in one space dimension for the nonlinear diffusion equation 

(1) dtu = in T'^, f > 0, m(0) = 

with negative exponent a < 0, where T'^ is the d-dimensional torus. Equation ([T]) can 
be written as the gradient flow of the entropy S[u] = {a{a — 1))“^ with respect 

to the Wasserstein distance. The gradient flow formulation gives rise to a natural semi¬ 
discretization in time of the evolution by means of the minimizing movement scheme, 
leading to a minimization problem for the sum of kinetic and potential energies. 

Minimizing movement schemes for evolution equations with an underlying gradient flow 
structure were first suggested by De Giorgi [12] in an abstract framework. Jordan, Kinder- 
lehrer, and Otto im have shown that the solution to the linear Fokker-Planck equation can 
be obtained by minimizing the logarithmic entropy in the Wasserstein space. Since then, 
many nonlinear evolution equations have been shown to constitute Wasserstein gradient 
flows, for instance the porous-medium equation [23], the Keller-Segel model [1], equations 
for interacting gases ra, and a nonlinear fourth-order equation for quantum fluids [14] . 

The minimizing movement scheme can be interpreted as an implicit Euler semi-discreti¬ 
zation with respect to the Wasserstein gradient flow structure. While for some time this has 
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been mainly used as an analytical tool, more recently numerical approximations of evolu¬ 
tion equations using this scheme have been proposed. In one space dimension, the optimal 
transport metric becomes flat when re-parametrized by means of inverse cumulative func¬ 
tions, which simplifies the numerical solution; see e.g. [Illll[2nil22]. For multi-dimensional 
situations, one approach is based on the Eulerian representation of the discrete solution 
on a fixed grid. The resulting problem can be solved by using interior point methods [8], 
hnite elements [7], or hnite volumes [9]. Another approach employs the Lagrangian rep¬ 
resentation, which is well adapted to optimal transport. Examples are explicit marching 
schemes [15], moving meshes [6], linear hnite elements for a fourth-order equation [13], re¬ 
formulations in terms of evolving diffeomorphisms nn, and entropic smoothing using the 
Kullback-Leibler divergence [21]. The connection between Lagrangian schemes and the 
gradient how structure was investigated in [TU]. In this paper, we will use the Lagrangian 
viewpoint. 

The minimizing movement scheme of De Giorgi is of hrst order in time only since it is 
based on the implicit Euler method. Concerning higher-order schemes, we are only aware 
of the paper [31]. There, second-order gradient how schemes were suggested for the Euler 
equations, with hnite diherences in space and the two-step BDF (Backward Diherentiation 
Formula) method or diagonally implicit Runge-Kutta (DIRK) schemes in time. 

In this paper, we propose a fully discrete second-order minimizing movement scheme 
using quadratic hnite elements in space and the two-step BDF method in time. We con¬ 
sider periodic point-symmetric solutions. The hnite-dimensional minimization problem, 
constrained by the mass conservation, is solved by the method of Lagrange multipliers 
which leads to a sequential quadratic programming problem. 

By construction, our numerical scheme is of second order both in time and space, it 
conserves the mass and dissipates the G-norm of the approximation of the quadratic 
ansatz functions at time step k, 

(2) ||(g‘+‘, g‘)||| = |(g*«)’"Af„g‘+‘ - 2(g‘+‘)Df„g‘ + l(g‘)D4g‘, 

where the matrix is dehned in the approximation of the Wasserstein metric on the 
space of the quadratic ansatz functions. We refer to Section [3] for details. It turns out that 
numerically, the relative G-norm decays exponentially fast to zero. Although we cannot 
expect for the multistep scheme that the discrete entropy decays exponentially fast, this 
holds true for the numerical experiments performed in this paper. Furthermore, also the 
discrete variance of the original variable u and the Lagrangian variable decay exponentially 
fast. 

The numerical tests also indicate that the decay rate of the entropy is increasing in the 
number of grid points, i.e., the discrete decay rates are smaller than the (expected) value 
for the continuous equation. This result is in accordance with the Endings of [21] for a 
finite-volume approximation of a one-dimensional linear Fokker-Planck equation. 

Let us briefly review the literature for the diffusion equation ([ID with a < 0. Equation ([ID 
with a = —1 appears in the modeling of heat conduction in solid Helium, where the solution 
u corresponds to the inverse temperature [26]. When this equation is considered in the 
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whole space or in a bounded domain with homogeneous Dirichlet boundary conditions, it 
is sometimes called the super-fast diffusion equation [2H1 Chap. 9]. The critical exponent 
of this equation in one space dimension is a = — 1. For a > — 1 and vP G we have a 

smoothing property, namely u{t) G for any f > 0 |2Hl Section 9.1]. For a < —1, no 

solutions exist with data in L^(M). The non-existence range in dimensions d > 2 contains 
even all negative exponents, a < 0 [2^. However, if d = 1 and a < — 1, there is a weak 
smoothing effect. Indeed, given G for p = (1 — a)/2, the solution u exists and 

is locally bounded in M x (0, cxd). Furthermore, if G Lp(M) then there is instantaneous 
extinction, i.e. u{t) = 0 in M for alH > 0 [28l Theorem 9.3]. 

Clearly, such results cannot be expected when the super-fast diffusion equation is con¬ 
sidered on the torus. We expect that global-in-time weak solutions exist, which converge 
to the constant steady state as t —)■ oo. Since we could not hnd any results on the existence 
and large-time asymptotics in the literature in that situation and since the use of negative 
exponents is less standard, we provide a (short) proof for completeness in the appendix; 
see Section 12.11 for details. 

The paper is organized as follows. The global existence result and the exponential 
decay of the solutions to the constant steady state as well as some basic properties of the 
Wasserstein distance for periodic functions are stated in Section |2l Section [3] is devoted to 
the description of the numerical scheme, and some numerical experiments are presented in 
Section 01 We conclude in Section 01 The appendix contains the proofs of the existence 
and large-time asymptotics theorems and the calculations of the coefficients of the matrix 
My, and the Hessian of the discrete entropy. 

2. Prerequisites 

2.1. Existence of solutions and large-time asymptotics. Equation ([T]) on the torus 
does not possess the non-existence or instantaneous extinction properties of the super¬ 
fast diffusion equation in the whole space since mass cannot get lost. In fact, we expect that 
for any a < 0, there exists a global weak solution. If the initial datum is nonnegative 
only, dl]) is still a singular diffusion equation. However, because of the fast diffusion, the 
solution becomes positive for all positive times and, by parabolic regularity theory, also 
smooth. 

Theorem 1 (Existence of weak solutions). Let a < 0 and let vP G L°°(T'^) satisfy > 0. 
If a = —1, we assume additionally that J^^^ogu^dx > —oo. Then there exists a unique 
weak solution to dH) satisfying u°' G L^(0,T; FTHT'^)), dtu G L^(0, T; ifhT'^)') for allT > 0, 
and 0 < m(x, t) < sup^d for x , t >Q. 

The proof of this theorem is based on a standard regularization procedure but we need to 
distinguish carefully the cases —1 < a < 0, a = —1, and a < —1. We present the (short) 
proof for completeness in Appendix For t ^ oo, the (smooth) solution converges to 
the constant steady state. Since this constant is positive, diffusion slows down when time 
increases. Therefore, we cannot expect instantaneous extinction phenomena. Still, we are 
able to prove that the convergence is exponentially fast with respect to the L^-norm. We 
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introduce the entropy 


Hs\u] = 


u 


^ dx — 


Ijd 


ijd 


udx] , (3 > 1. 


The steady state of ([T]) is given by u^o = fjd dx if vol(T'^) = 1. 

Theorem 2 (Exponential decay). Let u be a smooth positive solution to ([T]) and let 
vol(T‘^) = 1. Then 

|| m (^) -UooWmTd) < 

where Cs > 0 and for a < 0, 1 < B < 2, 

and Cb > ^ is the constant in the Beckner inequality (UnD; for -l<a<t), 13 = 2(1 - a), 

2(1-2a) 


A = 




In the hrst result, the decay rate A depends on sup^pd which seems to be not optimal. 
The decay rate in the second result depends on the L^-norm of only but we need 
a particular value of (3. The proof is based on the entropy method; see Appendix 
Stronger decay results have been derived for the fast-diffusion equation in the whole space 
or in bounded domains; see, e.g., [51128]. However, our proof is very elementary and just 
an illustration for the qualitative behavior of the solutions to ([T]). 


2.2. The Wasserstein distance for periodic functions. We recall some basic facts 
about mass transportation and Wasserstein distances following [13]. For more details, we 
refer to [a [30]. Let Af be a Riemannian manifold with distance d : X x X ^ M+ := [0, cxo) 
and let Pm{X) be the convex set of measures with hxed mass M > 0 on Af. The L^- 
Wasserstein distance of two measures yUi, /i 2 G Pm{X) is dehned by 


(3) 


hf"[hi,h2]^= inf 

7ren(/ii,/x2), 


d{x,yfdTr{x,y), 


’XxX 


where n(/ii,yU 2 ) is the set of transport plans connecting yi with /i 2 , i.e. the set of all 
measures on A x A with respective marginals yi and /i 2 , 

7r(A X A) = /ii(A), 7r{X X B) = y2{B) 


for all measurable sets A, B C X. If the measures yi and /i 2 possess densities ui and 
U 2 , respectively, with respect to a hxed measure on A, then we write, slightly abusing the 
notation, W[ui,U 2 \ instead of W[yi,y 2 \- 

In one space dimension, there exists an explicit formula to compute W. Let A = (a, h) C 
M be a (possibly inhnite) interval, and let /ii, /i 2 G Pm{X) be two measures. We dehne 
their distribution functions 


Ui : {a,b) ^ [0,M], Ui{x) = yi{{a, x]), i = l,2. 
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As these functions are right-continuous and monotonically increasing, they possess right- 
continuous increasing pseudo-inverse functions Gi : [0, M] —>■ [a, b], given by 

Gi{u) = inf{x e {a,b) : Ui{x) > u}, i = 1,2. 

Then [29], 

pM 

(4) W[ui,U2]^= / iGiiu)-G2iu)fdu. 

Jo 

This formula does not extend to X = T ~ (0,1) because of the topology induced by the 
periodic boundary conditions, d{x, y) = min{|a: — y\,l — \x — y\}. The reason is that mass 
can be transported either clock- or counter-clockwise (see [13] for details). However, if the 
densities Ui and U 2 are point-symmetric, (|1]) still holds. More precisely, let Ui{x) = Ui{l—x) 
for X G (0,1) and i = 1,2. Then (|3]) holds, where Gi : [0, M] —)■ [0,1] is the inverse function 
of Ui{x) = Ui{y)dy [TS] Lemma 2.2]. 

3. Time discretization and Lagrangian coordinates 

3.1. The semi-discrete BDF scheme. We introduce the second-order minimizing move¬ 
ment scheme. First, we explain the underlying idea for the hnite-dimensional gradient flow 

(5) X = —'V(j){x), t > 0, a;(0) = xq, 

where 0 : ^ M is a smooth potential. This equation can be approximated by the 

following minimization problem: 

= argmin<F(x), <F(x) = —||x — -1- 0(a:), 

where r > 0 is the time step size and x” is an approximation of x(nr). The minimizer 
x"^^^ is a critical point and thus, 

0 = = \x^+^ - x^) + V0(x”+^), 

r 

which corresponds to the implicit Euler scheme. 

Instead of the Euler scheme, we wish to discretize ([I]) by a multistep method. As an 
example, consider the two-step BDF (or BDF-2) method, 

1 0x^+2 - = -V0(x"+2)^ 

where x"" and are given. Writing this scheme as 

^(x" - x’^+2) - 2(x”+^ - x^+2) _ _^v0(x"+^), 

we see that is a critical point of the functional 

<h(x) = ——Ik"' — x|p + -|lx"^^ — xk + 0(3^)- 
4r r 
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More generally, the BDF-/c approximation of (j5|). 

k 

i=0 

for given x '^,can be formnlated as 

k—1 k—1 

-rV0(a:"+^) = ^ ^ ai{x^+^ - 0 :”+^), 

i=0 2=0 

since Yli=o ~ 0; minimization problem 

^ k—1 

(6) = argmin<f)(a;), $(a;) =- ai\\x'^~^^ — xW"^ + (j){x), 

2r ^ 

In a similar way, we may formnlate general mnltistep methods. We recall that BDF-/c 
schemes are consistent if Yli=o = 0; X]i=i = 1 ^^^1 zero-stable if and only if fc < 5 
[25| Section 11.5]. 

The same idea as above is applicable to gradient flows in the L^-Wasserstein distance. 
For this, we replace the norm by the Wasserstein distance. For eqnation ([1]), scheme 
([6]) tnrns into 

^ k—1 

(7) = argmin<F(M), <F(m) = ——''^aiW[u'^'^\uf + S[u], 

where S[u\ = {a{a — 1))“^ fj n" dx. Scheme ([7]) can be interpreted as a BDF-/c minimizing 
movement scheme. 


3.2. Lagrangian coordinates. Before introdncing the spatial discretization, we rewrite 
the scheme ([7]) in a more explicit manner, nsing the inverse distribntion fnnctions G and G* 
of u and u*, respectively, which were introdnced in Section ^I72\ The nnmerical procednre is 
similar to that in [T3] bnt onr higher-order scheme introdnces some changes. We call u = 
U(x) E [0, M] the Lagrangian coordinate, which is conjngate to the Enlerian coordinate 
x G T, and we refer to the inverse distribntion fnnction G as the associated Lagrangian 
map. For a consistent change of variables, we need to express the entropy S[u] in terms of 
the Lagrangian coordinates. With the formnla for the differential of an inverse fnnction, 

u{x) = d^U{x) = ^ , 

d^G{u) 

and the change of nnknowns x = G{u) nnder the integral in S[u], we obtain 


S[u\ = 


u{x) dx 


a{a — 1) Jj u{xy " a{a — 1) 


g{uj) ^du, 
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where g{oo) = Note that the exponent in the integrand is positive since a < 0. In 

terms of ( 7 , the expression for the Wasserstein distance in (jl]) becomes (see [121 Section 2.3]) 

/ M pM / pUJ \ 2 

{G{u)-G*{u)fdu = j yj ig{v) - 9*iv))dT]j du 


-M rM 


(M - max{r],r]'}){g{T]) - g*{v))igiv') 


g*{r]'))dr]dr]'. 


Jo Jo 

This expression is simply a quadratic form in g — g* and thus easy to implement in the 
numerical scheme. We summarize our results which slightly generalize Lemma 2.3 in [T3] . 


Lemma 3. Let the initial datum : T ^ M 6 e point-symmetric. Then the solution 
to the BDF-k scheme ((7)) is in one-to-one correspondence to the solution g"^^^ obtained 
from the inductive scheme 


( 8 ) 


gU+k _ ai-gjnin \[/(^)^ 

9 


with initial condition g^ = 1/u o G and given g^,... ,g"'~^^ ^ obtained from a lower-order 
scheme. The argmin has to be taken over all measurable functions g : [0, M] —)• (0, 00 ) 
satisfying the mass constraint g{u) du = 1 , and the function T is given by 


r-M pM 


k-1 


(9) vl;(77) = -- 


'0 ^0 


(M - max{r7, r]'}) ^ ai{g{r]) - dg dg 


i=0 


-M 


g{u) du 


cr(l ci) Jq 

Moreover, $(«) = '^{g) and the functions and g"' are related by 


( 10 ) 


X^= g^{g) dg. 

9^[^) Jo 


In Lagrangian coordinates, the problem has become a minimization problem in \l/(( 7 ) 
which is the sum of a quadratic form and a convex functional, hence it is convex. In 
the special case a = —1, the second integral in ([9]) is quadratic too which simplihes 
the numerical discretization. Therefore, we will consider mainly numerical examples with 
a = — 1 in Section m 


3.3. Spatial discretization. We approximate the inhnite-dimensional variational prob¬ 
lem (IHl) by a hnite-dimensional one. Minimization in (jH]) is performed over the hnite- 
dimensional space of quadratic ansatz functions. This generalizes the approach in [T3] . 
where only linear ansatz functions were used. We dehne the ansatz space as follows. 

Let N E N and a mesh {xq, ..., x^} on [0,1] be given with xq = 0 and xat = 1. Using 
OB, we construct the mesh Qn = {wq, Ci;i..., cutv} of [0,M]. Then uq = 0, un = M, 
Ui < cuj+i, and un-i = M — ut (point-symmmetry), where i = l,...,A^ — 1. Since we wish 
to introduce quadratic ansatz functions, we add the grid points Uj^i /2 = {ojj+i + Uj)/2 for 
j=0,...,iV-l. 
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The basis functions : T —)■ M are defined by 




(^Ar(a;) 


id—LJj — X 

Ld o Ld 'j _]_ 

Ldj^l—Ld 

0 


Ldi—Ld 

Ldl 

Ld—LdN — l 
M —LdN — I 
0 


for u e 

for oj E [uj, ojj+i], 
otherwise, 

for u E [0, Wi], 
for CO E [cjat-i, M], 
otherwise. 


j = 1,...,N -1, 


This set of piecewise linear functions is supplemented by the following piecewise quadratic 
basis functions: 


(j)N+j{uj) 


1 - 


/ 2 Ld—[idj — i-\-Ldj) \ 2 
\ Ldj—Ldj — i ) 


for UJ E [oJj-i, cUj], 
otherwise, 




The ansatz space is the set of all positive, piecewise quadratic functions g : [0, M] —)■ M+ 
of the form 

2N 

(11) g{u) = '^gj4>j{uj). 

i=i 


We call g := {gi, ..., g 2 N) G [0, oo)^-^ the associated weight vector. By dehnition of c()j, we 
have g{ujj) = gj for j = 1,..., N. Moreover, as g is point-symmetric, go = gN- 

Now, for given mass M > 0 and grid Hat C [0, M], we dehne the set C as the 

set of weight vectors g for which the associated interpolation g from flTTl) satishes the mass 
constraint. 


(12) 1 = ^ g{uj) duJ = J2 + \dN+}j - c^i-i). 

The Wasserstein metric for functions approximated in this way becomes 


-M rM 


2N 


2N 


W[u,u*Y = 


>0 Jo 

N 


{M - max{g,T]'}) '^{gi - g*)(j)i{v) ^{9j - 9j)(pj{v') dgdg' 

i=l j=l 

N 

'y ^ ^9i ~ 9i ){9j ~ 9j)(^ij + y ^ {9N+i ~ 9N+i){9j ~ 9j)bij 
=1 i,j=l 

N N 

y ^ {9i ~ 9i ){9N+j ~ 9N+j)^ji T y ^ (S'N+i ~ 9N+i){9N+j ~ 9N+j)'^iji 


*j'=l 


where 


*j'=i 


aij — 


*j'=i 


-M rM 


>0 Jo 


{M — max{? 7 , ?7'})0j(r7)0j(?7') dgdg', 
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I'M pM 

(13) bij = / / {M - max{r],r]'})(j)N+i{ri)(l>j{v') dvdr]', 

Jo Jo 

pM pM 

Cij= / {M-max{r],r]'})4>N+i{v)4>N+j{v')d'r]dr]'■ 

Jo Jo 

The coefficients a^, bij, and Cij can be computed explicitly. The explicit expressions are 
given in Appendix IB. II Setting A = {aij), B = (%), C = {cij), and dehning the matrix 

(14) = (M,,) = , 

we can formulate the above sum as 

2N 

win, u-f = W(9. - ».•)(» - s;). 

ti=i 

As the matrices A and C are symmetric, is symmetric, too. The matrix A corresponds 
to the linear approximation considered in [T3] . 


3.4. Minimization. The numerical scheme consists of the following hnite-dimensional 
variational problem; 

(15) = argmin TAr(g), 

SSGm 

k-1 2N 

where 4'iv(g) = -— ^ Mij{gi - g'^^^){gj - g]^^) + 

1=0 i,j=l 

and where g = (^fi,..., g 2 N), 


N-l 


( 16 ) 


■Swlgl = 


^^i+1 


a[a 


+ gN+i(t>N+i) doj 


1 ) ^ 

' i=0 


The functions TAr(g) and T( 5 f) from ([9]) are related by TAr(g) = T( 5 f) with a piecewise 
quadratic function g dehned from g by (fTTD . Since 4/^1 is convex for a < 0 and the set 
is convex, there exists a unique minimizer of fll5p . 


3.5. Fully discrete Euler-Lagrange equations. The minimizer gn+k in flT^ is subject 
to the mass constraint flT^ . by dehnition of the set G^. Therefore, instead of working on 
the set Gj^, it is more convenient to consider f[T5D as a constrained minimization problem 
for g on the larger set which is solved by the method of Lagrange multipliers A using 
the Lagrange functional 


L(g, A) = ^N(g) - 1 - 

i=i 


N 


ffj-1 + 5'i 2 . , . 

- 2 -^ ) (cjj - coj-i) 


2 
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A critical point of L satisfies the 2N conditions 


1=0 


2N 

i=l 


9?^') + 


dS 


N 


% 


j = l,...,27V. 


The precise valnes for dSj^/dgj are given in Appendix IB.21 for a = —1. The condition for 
the constraint is recovered from 


0 = G 


2Af+i 


dL 


N 

i-E 

i=i 


9j-i + 9j , 1 / N 

- 2 -^ gfi'N+i ) {^j 


The vector G[g, A] = (Gi,..., G 2 V+ 1 )''' ^ is the gradient of L(g, A) with respect 

to (g, A). We approximate a critical point nnmerically by applying the Newton method 
to the hrst-order optimality condition G[g, A] = 0. This leads to a seqnential qnadratic 
programming method, since at every Newton iteration step a qnadratic snbproblem has to 
be solved. 


3.6. Implementation. Let the solntion g at the nth time step be given and let g^°^ := g, 
A^°^ := 0. The iteration is as follows: 

g(^+i) gb) + (5g)(^+i)^ := A^*) + ((5A)("+^\ 

where is the solntion to the linear system 

Lf[g(*),A(*)]((5g)(*+'),(5A)(^+iy = -G[g(*),A«], 

where Tr[g^^\ A^^^] denotes the Hessian of whose entries are given in Appendix IB.21 
for a = —1. The iteration is stopped if the norm of (((5g)*'®’''^\ ((5A)^^'''^^) is smaller than 
a certain threshold (see Section H] for details). In this case, we dehne ;= g(^+i) and 
^n+i ._ yb+i) _|_ x)th time step. For the BDF-/c scheme, the valnes g^,..., g^“^ 

are computed from a lower-order scheme. In the numerical section below, we employ the 
BDF-2 scheme only such that g^ is calculated by the implicit Euler method. 

Note that the constrained minimization problem is exactly mass conserving by construc¬ 
tion, but the Newton iteration introduces a small error which depends on the tolerance 
imposed in the Newton method. 

In each iteration step, we need to invert the dense matrix H which is the sum of 
and the Hessian of S^. This is not a numerical challenge in the one-dimensional case we 
consider but it may become critical in multi-dimensional discretizations on hne grids. For 
a = —1, however, and the Hessian of S'at are constant matrices which signihcantly 
simplihes the Newton scheme. 

3.7. Choice of the initial condition. In order to compute the initial condition in La- 
grangian coordinates, we need to make precise the values g^ of the vector g G and 
the points of the spatial lattice which is moving as the solution evolves. Let the mesh 
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{xq, ... be given and set g^iojj) = l/u^{xj). Approximating the initial function by a 
linear ansatz function, we obtain 

(17) + ^°), j = l,...,N. 

i=l 

This is a system of linear equations in oii,... ,un. Choosing the uniform grid = j/N, 
we can solve this system explicitly. Indeed, since 

1 + g^), 

which can be solved for uj+ii 

Uj+i = Uj + ^{g^+1 + g^)~\ j = o,...,N -1. 

As is assumed to be point-symmetric, this is true for (oij) too. Finally, we approximate 
g^{u)j^i/ 2 ) by the arithmetic mean + g^), j = I,... ,N. Consequently, the weights 

gN+i vanish for i = 1,..., A^ in the expansion g^{u) = J 2 ‘^=o 9 j^ 3 y which is consistent with 
our approximation fITT)) . We also refer to the discussion in [T31 Section 2.8]. 

4. Numerical experiments 

In this section, we present some numerical results for ([T]) with a = — 1, by employing the 
BDF-2 method with quadratic ansatz functions. We choose a uniform grid for x G [0,1] 
with N = 100 grid points, and the time step size r = 10“^. The Newton iterations are 
stopped if both the relative i°° error in the gf-variables and the norm of are 

smaller than 10“®. 

Figure [T] illustrates the temporal evolution of the solution u{x,t) with the initial condi¬ 
tions u^{x) = cos(27rx)^ -|- 0.01 (left hgure) and u^{x) = ^\x — 0.5| -I- 0.0001 — 0.1 (right 
hgure). We observe that, as expected, the solutions converge to the constant steady state. 
Because of the negative exponent a, very small initial values increase quickly in time. 

A nice feature of the Wasserstein gradient flow scheme is that we may interpret the 
evolution as a process of redistribution of particles with spatio-temporal density u{x,t) on 
T under the influence of a nonlinear particle interaction, which is described by S. The 
way in which the initial density is “deformed” during the time evolution is illustrated 
in Figure [2l We have chosen 50 “test particles” for the solutions to ([T]) for the initial 
conditions chosen above. We stress the fact that the density of trajectories can generally 
be not identihed with the density u of the solution. 

We verify that the discretization is indeed of second order. Figure [3] shows the i°°- 
error for various numbers of grid points N. We have chosen the initial datum Uq{x) = 
cos(27rx)^ -|- 0.1, the end time T = 0.004, and the time step size r = 10“^. The reference 
solution is computed by using N = 500, and r = 10“^. The differences g{-,T) — g^ei^-^T) 
and u{-,T) — Urei{-,T) in the £°° norm feature the expected second-order dependence on 
N. 
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X-space 


Figure 1. Time evolution of the solution to the diffusion equation ([T]) with 
a = —1 for two different initial conditions. 




X-space X-Space 

Figure 2. Particle trajectories in the Wasserstein gradient flow scheme, 
corresponding to the solutions of Figure [T] with iV = 50. 


Next, we hx the number of grid points N = 100 and compute the L^(T)) error 

for varying time step sizes r; see Figure IH Because of the approximation of the initial 
datum as detailed in Section 1X71 the error will be not of second order initially. Therefore, 
we compute the error in the interval (r*,T) with r* = 10“^. The errors are of second 
order, as expected. 

The time decay of the discrete version of the relative entropy S[u\ — S'[moo] is presented 
in Figure |5] (left) for various grid numbers. We observe that the decay is exponential until 
saturation. The saturation comes from the spatial error and the error from the Newton 
iteration. The decay rate is estimated in the linear regime from the difference quotient 

A ^(logS'[M(t + r)] - logS'[M(t)]). 
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Number of grid points 


Figure 3. £°°-error for {g — g-ref)(T) (left) and {u — u^e{)(T) (right) at T = 
0.004 for various numbers of grid points. 




Time step size Time step size 

Figure 4. f'°°(r*,T;£^(T))-error for g — g^-ef (left) and u — Uref (right) for 
various time step sizes, with t* = 10“^. 

The numerical decay rates for a = —1, —2 are shown in Figure [6l The rate for a = —2 
(right) is much larger than the corresponding one for a = —1 (left), since a smaller exponent 
yields a larger diffusion coefficient (if m < 1) and thus, diffusion becomes faster. We also 
see that the decay rates become larger on a hner spatial grid. This behavior seems to 
conhrm recent analytical results for spatial discretizations of Fokker-Planck equations; see 
pn Section 5]. One may ask if a similar behavior can be observed for the decay rate 
as a function of the time step size. However, our numerical experiments do not show a 
monotonic dependence (hgures not presented); rather the decay rates vary in a small range 
which seems to be determined by the other numerical error parts. 

In Figure [71 the decay of the square of the relative G-norm is presented. The G-norm is 
calculated according to ([2]), where the argument is given by g — goo and goo is the weight 
vector corresponding to the constant steady state. Again, the decay is much faster for 
a = —2 because of the faster diffusion. 
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Figure 5. Discrete relative entropy S'[M(t)] — S'[moo] versus time with a = —1 
for various N. 






Number of grid points Number of grid points 


Figure 6. Estimated decay rate of the entropy versus number of grid points 
TV for a = —1 (left) and a = —2 (right). 


Finally, we present some results on the time decay of the discrete variance of and 
at time rn, dehned by 


N+l 

Var(M”)^ = - Ef{xi - Xi-i), 

1 = 1 

where E is the expectation value of (which equals the mass and is therefore constant in 
time). The discrete variance of is dehned in a similar way. Interestingly, the variances 
are exponentially decaying (Figure |8]), although it is not clear how to prove this property 
analytically. 
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Figure 7. Relative G-norm versus time for various grid numbers N (left) 
and for two values of the exponent a (right). 



Figure 8. Discrete variance of n” (left) and g"' (right) versus time for var¬ 
ious grid numbers N. 


5. Conclusion 

We have proposed an extension of the standard minimizing movement scheme to higher- 
order BDF time discretizations. Quadratic hnite elements have been used to discretise in 
(mass) space. As an example, we have considered a (singular) diffusion equation, but it 
should be possible to adapt the scheme to other partial differential equations which consti¬ 
tute L?‘ Wasserstein gradient flows. The implicit numerical scheme is based on successive 
solution of constrained minimization problems. The quadratic convergence in space and 
time of our method has been conhrmed numerically. It turns out that the relative entropy, 
G-norm, and variance converge exponentially to zero. An interesting observation is that 
the decay rate depends monotonically on the number of grid numbers, but it seems that 
there is no monotonic behavior with respect to the time step size. Future work may be 
concerned with an analytical derivation of numerical decay rates using a discrete variant 
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of the Bakry-Emery approach. First steps in this direction were presented in [21] bnt only 
for a semi-discretization. Possibly this approach has to be adapted to G-stable schemes in 
the spirit of [T8] . 


Appendix A. Proofs of Theorems □ and [2] 


Proof of Theorem [H The existence proof is based on a regnlarization of the initial datnm 
and a hxed-point argnment. Let 0 < £ < 1 and Mo,e = + s. Let Qt = x (0,T) and 

M = snp-jpd Set 

K = {u E L‘^{Qt) : e < U < M, ||M||L 2 ( 0 ,T;Rl(T'i)) + ||^i'*^||L 2 ( 0 ,r;Hl(T'i)') < C}, 

where C > 0 will be determined later. The set K is convex and, by Anbin’s lemma, 
compact in Let v E K and let u E L^{Qt) be the weak solntion to 

(18) dtu = div(n“-VM) in t > 0, m(0) = 

This dehnes the hxed-point operator Z : K ^ L'^{Qt), v ^ u. Standard argnments 
show that Z is continnous. We verify that Z{K) C K. By the maximnm principle, 
e < u < M. Using n as a test fnnction in the weak formnlation of (fT8|) shows that 
ll'^llL2(o,r;Hi(T‘i)) ^ C'i(e), where Ci{e) > 0 is some constant depending on e. Moreover, 
ll^tw||L 2 ( 0 ,T;rii(T‘*)') < < C 2 {e). Thns, setting, C := Ci{e) + C 2 {e), we infer 

that u E K. By the hxed-point theorem of Schander, there exists a hxed point Us of Z. 

In order to perform the limit £ —?■ 0, we need to derive ^-independent estimates for Ue- 
To this end, we need to distingnish several cases. First, let a = —1. Employing the test 
fnnction 1 — 1/u^ in the weak formulation of (ITSD with u = v = u^, we hnd that 


Ijd 


{Us{t) - log Me(t)) dx + 



0 


|VMg f dxds = / -f £ — log(M^ -|-£)) dx. 


Ijd 


The right-hand side is uniformly bounded as we assumed that — logw^dx < oo. Since 

we infer uniform estimates for Ue in L^(0, T; iL^(T‘^)) and also in 

Next, let a 7 ^ —1. The test function uf in the weak formulation of 0181) gives 


o -)-1 


dx H— 



ifd a Jq jfd 

If —1 < a < 0, we write this equation as 


I dxds = 


ijd 



0 


If a < —1, we obtain 
1 


I Vm"|^ dxds = 


a 


< 


d T 1 jTd. 

a f 


o -|- 1 


W(i) dx + 
dx. 


o -|- 1 


a 


U T 1 dfi 


dx. 


{u^ -F dx 


Ij’d 


f Me(t)“+^dx + — r [ 

T'i Jq Jjd 


Vm“|^ dxds = 


-a 


Ijd 


{u^ + e)“+^ dx. 


—a — 1 
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In both cases, since is assumed to be bounded, we infer a uniform bound for 
in L^(0, T; and consequently also for dtU^ in L^(0, T; Moreover, in 

view of |Vm“P = |Vwep, it follows that (ug) is bounded 

in L2(0,T; 

We infer for all a < 0 the following bounds: 

l|w|lL°°(0,r;L°“(T‘i)) + ll'^e |lL2(o,r;Hi(T‘i)) + \We\\L‘2{0,T-,H^{T<i)) + \\dtUe\\L'^{0,T-,H^{T‘^y) < C 3 . 

By Aubin’s lemma, there exists a subsequence which is not relabeled such that ^ u 
strongly in LF‘{Qt) as £ —)■ 0. Moreover, m" ^ weakly in L^(0, T; and dtUf. 

dfU weakly in L^(0, T; hf^(T'^)'). Thus, we may pass to the limit e —)■ 0 in the weak 
formulation which shows that u solves ([1]). □ 


Proof of Theorem [H Employing ([T]) and integration by parts, we hnd that 


dH, 




dt 


= -P{P-l) w“+/3-3|Vn|2c^a: =-|(/3-l) / dx. 


ijd jd jfd 

We employ the bound u < M = sup^d and the Beckner inequality [3] (note that we 
assumed that vol(T°') = 1), 

h 

(19) I 

JT<i 

to obtain 

dHp 

dt 


^dx - ( f udx^ <Cb [ \Vu^/^\^dx for G H\n), 1 < /3 < 2, 

4(/3-l) 


fdCB 


^ dx 


Ijd 


u dx 


ijd 


fdCB 




Then Gronwall’s lemma yields H[u{t)] < H[vf‘]e with A = —4(/5 — 1)M" ^/(/3 Gb) for 
l<(3<2. 

For the second result, let —1 < a < 0 and fd = 2(1 — a). Similarly as above, we hnd that 
dHfi _ 4/3 (/3 - 1) 


dt 


< - 


{a + (d - 1)2 
8(1-2a) f 


ijd 


|y^(a+/3-l)/2|2^^ 


1 — a 
8(1-2a) 
1 — a 


u 


dx - 


JTd. \ Jjd 

H^/2[u\, 


p \ a+/3-l 

/ udx ] 
hd J 


since 
\\u 


a + /3 — 1 = /3/2 G (1,2]. Using the inequalities ||M^/^||i 2 (']rd) > ||M^/^||ii(Td) and 
= l|■*^llL]8^(Td) — (again we employ vol(T'^) = 1 here), it follows that 

-^/3M = IIl2(T'') “ II'*^IIli(t'^) 

= (ll“^^^llL2(Td) + \\u\fLi(jd-^ {\\u^^‘^\\L'^(Td) - ll'wlliqTd)) 


> {Wu^^’^Wh^Td) + ||M||^{(Td))(||M^^lLi(T<^) - \\u\\%\jd)) 


1/3/2 


/3/2| 


|/3/2 
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,/3/2|| II,, 11/3/2 ^ r,,l ^ oIL,I|/3/2 


Since the solution to ([T]) conserves mass, ||M(t)||ii(Td) = ||m°|| and we end up with 

dHp 


4(1 -2a) 

dt - n _ ..mu.oii/3/2 


< 


(1 -«)l|w°llLi(Td 


-HM 

) 


Then Gronwall’s lemma shows that Hi 3 [u{t)] < Hp[u^]e with A = 4(1 — 2a)/((l — 


a)\\u 


On 1—0 


Finally, the statement of the theorem follows after applying the generalized Csiszar- 
Kullback inequality in the form 


\u — V 


< Cf}\\v\\L^jd) ^udx 


1 < /3 < 2. 


for functions u, v E L^(T‘^) such that Jr^^udx = Jr^^^dx. The proof is a slight generaliza¬ 
tion of the proof of Theorem 1.4 in [16] taking ip(t) = □ 

Appendix B. Computations 

In this appendix, we detail the calculations for the coefficients of the matrix flT^ . and 
the Hessian of the discrete entropy flT6|) . both in the case a = — 1. 

B.l. Computation of the coefficients Mij. We compute the coefficients of the matrix 
JH, u. the coefficients Oij, bij, and Cij dehned in flT^ . In the following, we set 

~ = 2 ~ ^i-i)’ = 2 (^i+i + 

Lemma 4 (Coefficients aij). The coefficients of the symmetric matrix A = {aij), defined 
in ca, read as 

ajj = A]{M - aj) - ^(12A] + 6] + 1 < j < iV - 1, 


tj+i ~ ^j+i) 


5? 


1 < J < A - 1, 


oj+i ~ 220’ 

Ojk = AjAk{M - (Tfc), j + 2 < k < N - 1, 

A3 

’ 120 ’ 


1 f UJo 

aiN = -AiAat \ M —— 

= 2^/^/v (^M - cTj , 


2<j<N-2, 

Al 


'^N-l,N — 2^/V-lAAr^M — -{un-2 + 2un-i) ) — 

a - ^a2 I 

^NN - 


N 


120 ’ 
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This lemma was proved in [131 Lemma 2.6]. For the convenience of the reader, we recall 
the proof below. 


Proof. We reformulate the integral in flT^ : 


-M 


fM 


-M 


-M 


ajk = M / (l)j{ri)dr] / (fkiv') dr]'- / (l)j{r])dr] / T]'(j)kiv') dr]' - Jjk, 

Jo Jo Jo Jo 

pM pM 

where Jjk= / / {r] - r]')+(j)j{r])(j)k{r]') dr]dr]', 


>0 JO 


where [r] — 77 ')+ := max{0, r] — r]'}. The hrst two integrals become 


fM 


(j)j{r]) dr] = 6j, 


r]'(j)k{r]') dr]' = A^o-fc. 


Since A is symmetric, it is sufficient to consider 1 < j < k. If j + 2 < k < N, the support of 
4>j{r])(t>kij]') is contained in oij+i] x [uk-i x 01 ^+ 1 ]. Hence, the support is nonvanishing 
if 77 < oij+i < oJk-i < ? 7 ^ but then (77 — 77 ')+ = 0 except for 77 = 77 '. We conclude that 
Jjk = 0 and it is sufficient to compute only Jjj and Jjj+i- 




djj — 


dj,j+i — 


iiv) 


rr) 


' — 1 \ — 1 

07 (^) 


(77 - 77 ') 0 ,( 77 O dri ) dr] = 


lUj-l 




\ 6^ 

iv - v')(pj{v') dv'j = 


Next, let k = N. Then the support of (j)N is contained in [0,a;i] U [u)n-i,M], and we 
compute: 


fM 


fM 


-M 




OjN = M / (t)j{j])dr] / 4)N{r]')dr]' 


r](j)j{r]) dr] / cjrNir]') dr]' 


fM 


fM 


(pj{r])dr] 


V'MV') dr]' - K+ - K. 


I UlN-l 


= M -a,+ 


’N 


- Kt - K: 


where 


77/:= 


K- := 


pM pull 

<0 A 

pM pM 


W - dr]dr]', 


'0 JuiN-l 


iv - V)+ 07 (^) 0 v(V) dr]dr]'. 


For 2 < j < N — 2, the supports of (pj and 0Ar do not intersect such that Kj" = 0. For 
i = 1, we have Kf = 0 and Kf = A^/120, whereas for j = A — 1 , = 0 and 
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= A|f/120. Furthermore, = A|^/30. Moreover, since 5^^ = ujn — <^v-i = 
(M — cuo) — {M — cui) = Wi, 

M-ai + ^= M- ^(cuo + wi + U2) + = M - y. 

Collecting these results, the lemma follows. □ 

Lemma 5 (Coefficients bij). The coefficients of the matrix B = {bij), defined in flT^ . read 
as 


^j^j) ffi 1 A J A -^5 

~ 2 ~ i^j+2 ~ ^ — j A ~ 1, 

bjk = ‘^{MdjAk - SjdkUk), I < j < k < N, 

bjk = ^ {2MdjAk - - cnj) A^), j > k - 2, 

2 1 1 

blN = — -r{^2 ~ ^l)^N — -^{‘2 uJn ~ ^"n-I ~ — filN■, 

00 9 

2 1 1 

bjN = — -(<^1+1 — ^^j)dN — '^{2,ulf — — u)N_iU)N)6j, 2 < j < N — 1, 


9 


3 

2 1 2 2 1 2 2 

bNN = -^nAn — -{^N ~ ^N-i)^N — ffi2,Uj^ — CJjv-l ~ ^N-l^N)dN — PnN, 
00 9 


where 


n ^2 ^2 ^2 ^ 2^ 


1 7 11 

g^i + 2^i + l + - g^i^i+2 + g^i^j + 2, 


/^i+i,i - ^ (^1+2 ~ ^i+i + ^i‘^'j+i^j+2 - a;|+2^i+i)) > 
fiiN = ~ 1^1 + 3 ( 022^1 — 

PnN = ^(l^AT+i — l^N + ^{^N^N+1 — (^n+1^n))- 

Proof. The computation is similar to the previous proof. We write the integral for bjk with 
j, A; < A — 1 and for j < k as 


-M 


-M 


r-M 


-M 


bjk = M / (l)N+j{r]) dr] / fikij]') H 

Jo Jo 

and for j > k as 


fiN+jijl) dr] / v'Mv') dr]' - Jj^, 


-M 


-M 


r-M 


-M 


bjk = M (pN+jiv)dv (pk{v')dv'- T](pN+jiv)dv fikiv) dv'- 
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where 


Jjk - 


j+ — 
^jk - 


pM pM 

lo Jo 

pM pM 

lo Jo 


iv - v')+(l>N+j{v)Mv') dr]dr]', 
W - v)+(I^N+j{v)Mv') dvdv'- 


We compute 


fM 


-M 


(j)N+j{v) dv = 


V(j)N+j{v) dr] = 


The integrals vanish ii k ^ j — 1, j since (r] — r]')^ and [j]' — 17 )+ vanish. This proves 
the expressions for with j < k — 1 and j > k + 2. The coefficients bjj and bj^ij are 
calculated in the same way. 

It remains to compute the matrix coefficients coming from the boundary elements. The 
computation of bjN is straightforward as the support of 02 v is contained on the single 
subinterval [un-i, M]. For the boundary elements bj^, we take into account that the 
support of (j)N is contained in [u)n-i,M] and [0,a;i], which yields 


"M 


i-M 


»M 


bjN = M (pN+jir])dr] (pNiv') dr]' 
Jo Jo 

where 


r](j)N+j{r]) dr] 4>N{r]') dr]'- - I^nn, 


fJlN — 


/Jnn — 


rM PLOi 

'0 Jo 

rM rM 


iv' - ^)+0N+j(h)0N(h') dr]dr]', 


iv - v')+(l>jiv)(l>N{v') dr]dr]', 


JO Jujrf-1 

and the computation is as before. 


□ 


Lemma 6 (Coefficients Cij). The coefficients of the symmetric matrix C = {cij), defined 
in TO, read as 

9*^^ 9 


o 2 2 

~ — (^j) — —{(^j+i — (^j — — ^j))i 1 < J ^ 


Cjk — gdjdk ^ ^ j < k < N. 


Proof. We compute 


fM 


-M 


-M 


-M 


Cjk = M (pN+jir])dr] (pN+kiv') dr]’ 
Jo Jo 

pM pM 

where ^jk = 


(pN+j{v)dv / v'4>N+k{v')dv' -Ijk, 


iv - v')+(l>N+jiv)(l>N+kiv') dvdi- 


'0 JO 
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As before, we find that 'jjk = 0 for all j 7 ^ k. Moreover, 


= / (pN+jiv) dv I iv- v')<pN+j{v') dr! = — - u] - 3{uj+iUj{uj+i - uj)) , 


JUJj-l J 

which hnishes the proof. 


□ 


B.2. Computation of the coefficients of the Hessian of Aat. We compute the gradient 
and Hessian of the discrete entropy flT^ for the case a = —1. We set for k = 0,..., N — 1: 

1 r‘^k+1 

>S'v,fc[g] — 7 ^ {9k(pk{^) + gk+l(pk+l{^) + gN+k(pN+k{^)) du, 

2 Jujf, 

where g = {gi,...,g 2 N) G Furthermore, we abbreviate dkSjsfj = dSNj/dgk and 

dj,kSN,£ = dSN^i/dgjdkg. A computation shows that 

dkSN,k[s\ = “^(2(5'v+fc + gk) + 5'fc+i)) 

2 

dkSN,k-i[s\ = -^dk{gk — gk-i — gN+k-i), 

= dk+i + gk) + • 

As SN,k and Sisf,k-i depend on gk, we obtain (recall (1T6|) ) 

dkSN = dkSN,k + dkSN,k-i, d^+kSN = dN+kSN,N+k- 
The second-order derivatives become 


- --4. 


9t.kS^.t-i - 34 . 


- i4+i, - - 54 . 

9N-\-k,k+l^N,k dN+kyN+k^Nyk ~^dk-\-\- 

3 15 


Then the elements of the Hessian of Sj^ read as 


dkykSNyk — 
dkyN+kSNyk = 


dkyk-l^N = dkyk-l^Nyk-l, 
dkyk+lSN = dkyk+lSNyk, 
dk,N+k-lSN = dkyN+k-lSNyk-1, 
dN + kyN + k^N = dN + kyN + kSN , k - 


dkykSN = dkykSNyk + dkykSNyk-l, 
dk,N+kSN = dkyN+kSNyk, 
dN+kyk+lSN = dN+kyk+lSNyk, 
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